1 Import library

Code

# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)

# table
library(gt)

# map
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle
library(oce)

# data wrangling
library(tidyr)
library(dplyr)
library(purrr)
library(forcats)
library(lubridate)

2 Setting up custom function

2.1 `windrose`

Code

windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet = F) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, 7, 20), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin lables
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet) {
      data2$origin <- unique(data_to_plot$origin)
    }

    # set up max value
    maxvalue <- 35

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "white",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code

# import the data
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartureLocation == "ANNU")

Code

# get solar angle when data location
solar_angle_inter <- data_dive %>%
  filter(!is.na(Lat)) %>%
  mutate(solar_angle = sunAngle(date, Long, Lat)$altitude) %>%
  select(id, DiveNumber, solar_angle)

# merge that with data_dive
data_dive <- merge(data_dive,
  solar_angle_inter,
  by = c("id", "DiveNumber"),
  all.x = T
)

# add day-night
# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/
data_dive <- data_dive %>%
  mutate(day_night = if_else(solar_angle < -18, "night", "day"))

Let’s add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, i.e. if this difference is above xxx m then is not considered as a benthic dive.

Code

data_dive <- data_dive %>%
  # change DiveType to include condition on difference
  # between max depth and bathymetry (<50m)
  mutate(
    DiveType_50 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 50, 0, DiveType),
    DiveTypeName_50 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 50,
      "Transit",
      DiveTypeName
    ),
    # same with <100m
    DiveType_100 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 100, 0, DiveType),
    DiveTypeName_100 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 100,
      "Transit",
      DiveTypeName
    ),
    # same with <150m
    DiveType_150 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 150, 0, DiveType),
    DiveTypeName_150 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 150,
      "Transit",
      DiveTypeName
    ),
    # same with <200m
    DiveType_200 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 200, 0, DiveType),
    DiveTypeName_200 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 200,
      "Transit",
      DiveTypeName
    )
  )

4 Maps

Code

# check if background_ggoceanmaps exist
if (file.exists("../export/background_ggoceanmap.rds")) {
  trip <- readRDS("../export/background_ggoceanmap.rds")
} else {
  # using ggOceanMaps
  trip <- basemap(
    limits = c(170, -110, 20, 59),
    bathymetry = TRUE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip <- trip +
    geom_sf(data = graticule, color = "grey50")

  trip$layers <- trip$layers[c(1, 3, 2)]

  # save result
  saveRDS(trip, "../export/background_ggoceanmap.rds")
}

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "All Benthic Dives",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

Figure 1: Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_50 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<50",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

Figure 2: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_100 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<100",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

Figure 3: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_150 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<150",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

Figure 4: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_200 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<200",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

Figure 5: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

5 Windrose

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(3000, 6000, 9000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 6: Circular histogram plots displaying the time when they perform benthic dives across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.5 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 7: Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_50 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 8: Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_100 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 9: Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_150 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 10: Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_200 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

6 Boxplot of the benthic dive rate (day vs. night)

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop") %>%
  group_by(id, day_night) %>%
  summarise(nb_average_benthic_dive = mean(nb_benthic_dive), .groups = "drop")


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of dives per individual per hour"
  )

Figure 11: Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time

Code

# Mann Whitney / Wilcoxon rank sum test on benthic dive rate
broom::tidy(wilcox.test(
  data_plot %>% filter(day_night == "day") %>% pull(nb_average_benthic_dive),
  data_plot %>% filter(day_night == "night") %>% pull(nb_average_benthic_dive)
)) %>%
  gt() %>%
  tab_header(title = "Mann Whitney / Wilcoxon rank sum test on benthic dives rate")

statistic	p.value	method	alternative
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
62027.5	8.80423e-09	Wilcoxon rank sum test with continuity correction	two.sided

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop") %>%
  group_by(id, day_night) %>%
  summarise(nb_average_benthic_dive = mean(nb_benthic_dive), .groups = "drop")


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of dives per individual per hour"
  )

Figure 12: Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time

Code

# Mann Whitney / Wilcoxon rank sum test on benthic dive rate
broom::tidy(wilcox.test(
  data_plot %>% filter(day_night == "day") %>% pull(nb_average_benthic_dive),
  data_plot %>% filter(day_night == "night") %>% pull(nb_average_benthic_dive)
)) %>%
  gt() %>%
  tab_header(title = "Mann Whitney / Wilcoxon rank sum test on benthic dives rate")

statistic	p.value	method	alternative
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
94857.5	9.608421e-07	Wilcoxon rank sum test with continuity correction	two.sided

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop") %>%
  group_by(id, day_night) %>%
  summarise(nb_average_benthic_dive = mean(nb_benthic_dive), .groups = "drop")


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of dives per individual per hour"
  )

Figure 13: Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time

Code

# Mann Whitney / Wilcoxon rank sum test on benthic dive rate
broom::tidy(wilcox.test(
  data_plot %>% filter(day_night == "day") %>% pull(nb_average_benthic_dive),
  data_plot %>% filter(day_night == "night") %>% pull(nb_average_benthic_dive)
)) %>%
  gt() %>%
  tab_header(title = "Mann Whitney / Wilcoxon rank sum test on benthic dives rate")

statistic	p.value	method	alternative
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
94412.5	1.916629e-06	Wilcoxon rank sum test with continuity correction	two.sided

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop") %>%
  group_by(id, day_night) %>%
  summarise(nb_average_benthic_dive = mean(nb_benthic_dive), .groups = "drop")


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of dives per individual per hour"
  )

Figure 14: Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time

Code

# Mann Whitney / Wilcoxon rank sum test on benthic dive rate
broom::tidy(wilcox.test(
  data_plot %>% filter(day_night == "day") %>% pull(nb_average_benthic_dive),
  data_plot %>% filter(day_night == "night") %>% pull(nb_average_benthic_dive)
)) %>%
  gt() %>%
  tab_header(title = "Mann Whitney / Wilcoxon rank sum test on benthic dives rate")

statistic	p.value	method	alternative
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
94853.5	9.670154e-07	Wilcoxon rank sum test with continuity correction	two.sided

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop") %>%
  group_by(id, day_night) %>%
  summarise(nb_average_benthic_dive = mean(nb_benthic_dive), .groups = "drop")


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of dives per individual per hour"
  )

Figure 15: Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time

Code

# Mann Whitney / Wilcoxon rank sum test on benthic dive rate
broom::tidy(wilcox.test(
  data_plot %>% filter(day_night == "day") %>% pull(nb_average_benthic_dive),
  data_plot %>% filter(day_night == "night") %>% pull(nb_average_benthic_dive)
)) %>%
  gt() %>%
  tab_header(title = "Mann Whitney / Wilcoxon rank sum test on benthic dives rate")

statistic	p.value	method	alternative
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
94946.5	1.710967e-06	Wilcoxon rank sum test with continuity correction	two.sided

7 Boxplot of the number of benthic dive (day vs. night)

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 1000)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 16: Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	5.96	0.05	117.16	0.00
day_nightnight	−0.84	0.09	−9.07	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 17: Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.54	0.07	65.08	0.00
day_nightnight	−0.54	0.12	−4.69	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 18: Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.61	0.07	65.55	0.00
day_nightnight	−0.54	0.12	−4.65	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 19: Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.64	0.07	65.52	0.00
day_nightnight	−0.54	0.12	−4.63	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 20: Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.67	0.07	65.47	0.00
day_nightnight	−0.55	0.12	−4.66	0.00

1 Import library

2 Setting up custom function

2.1 windrose

3 Import data

4 Maps

5 Windrose

6 Boxplot of the benthic dive rate (day vs. night)

7 Boxplot of the number of benthic dive (day vs. night)

2.1 `windrose`